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1.  Introduction 


1.1  Description  of  Problem 

Dynamic  failure  in  bound  particulate  materials  is  a  combination  of  physical  processes  including 
grain  and  matrix  deformation,  intra- granular  cracking,  matrix  cracking,  and 
inter-granular-matrix/binder  cracking/debonding,  and  is  influenced  by  global  initial  boundary 
value  problem  (IB VP)  conditions.  Discovering  how  these  processes  occur  by  experimental 
measurements  is  difficult  because  of  their  dynamic  nature  and  the  influence  of  global  boundary 
conditions  (BCs).  Global  BCs,  such  as  lateral  confinement  on  cylindrical  compression 
specimens,  can  influence  the  resulting  failure  mode,  generating  in  a  glass  ceramic  composite  axial 
splitting  and  fragmentation  when  there  is  no  confinement  and  shear  fractures  with  confinement 
(4).  Thus,  we  resort  to  physics-based  modeling  to  help  uncover  these  origins  dynamically. 

Examples  of  bound  particulate  materials  include,  but  are  not  limited  to,  the  following: 
polycrystalline  ceramics  (crystalline  grains  with  amorphous  grain  boundary  phases,  figure  1(a)), 
metal  matrix  composites  (metallic  grains  with  bulk  amorphous  metallic  binder,  figure  1(b)), 
particulate  energetic  materials  (explosive  crystalline  grains  with  polymeric  binder,  figure  1(c)), 
asphalt  pavement  (stone/rubber  aggregate  with  hardened  binder,  figure  1(d)),  mortar  (sand  grains 
with  cement  binder),  conventional  quasi-brittle  concrete  (stone  aggregate  with  cement  binder), 
and  sandstones  (sand  grains  with  clayey  binder).  Bound  particulate  materials  contain  grains 
(quasi-brittle  or  ductile)  bound  by  binder  material  oftentimes  called  the  “matrix.”  The 
heterogeneous  particulate  nature  of  these  materials  governs  their  mechanical  behavior  at  the 
grain-to-macro-scales,  especially  in  IB  VPs  for  which  localized  deformation  nucleates.  Thus, 
grain- scale  material  model  resolution  is  needed  in  regions  of  localized  deformation  nucleation 
(e.g.,  at  a  macro-crack  tip,  or  at  the  high  shear  strain  rate  interface  region  between  a  projectile  and 
target  material1).  To  predict  dynamic  failure  for  realistic  IB  VPs,  a  modeling  approach  will  need 
to  account  simultaneously  for  the  underlying  grain-scale  physics  and  macro-scale  continuum 
IB  VP  conditions. 

Direct  Numerical  Simulation  (DNS)  represents  directly  the  grain-scale  mechanical  behavior 
under  static  (5)  and  dynamic  loading  conditions  (6-8).  Currently,  DNS  is  the  best  approach  to 
understanding  fundamentally  dynamic  material  failure,  but  is  deficient  in  the  following  ways: 


'Both  projectile  and  target  material  could  be  modeled  with  such  grain-scale  material  model  resolution  at  their 
interface  region  where  significant  fracture  and  comminution  occurs. 
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Figure  1.  (a)  Microstructure  of  alumina,  composed  of  grains  bound  by  glassy  phase  (Sandia).  (b)  SiC  rein¬ 
forced  2080  aluminum  metal  matrix  composite  (i).  The  four  large  black  squares  are  indents  to  identify  the 
region,  (c)  Cracking  in  explosive  HMX  grains  and  at  grain-matrix  interfaces  (2).  (d)  Cracking  in  asphalt 
pavement. 
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(i)  it  is  limited  by  current  computing  power  (even  massively-parallel  computing)  to  a  small 
representative  volume  element  (RVE)  of  the  material;  and  (ii)  it  usually  must  assume  unrealistic 
BCs  on  the  RVE  (e.g.,  periodic,  or  prescribed  uniform  traction  or  displacement).  Thus, 
multi-scale  modeling  techniques  are  needed  to  predict  dynamic  failure  in  bound  particulate 
materials. 

Current  multi-scale  approaches  attempt  to  do  this  but  fall  short  by  one  or  more  of  the  following 
limitations:  (i)  not  providing  proper  BCs  on  the  micro-structural  DNS  region;  (ii)  homogenizing 
at  the  macro-scale  the  underlying  micro-structural  response  in  the  unit  cell  and,  thus,  not 
maintaining  a  computational  ‘open  window’  to  model  micro- structurally  dynamic  failure2;  and 
(iii)  not  making  these  methods  adaptive,  i.e.,  moving  a  computational  ‘open  window’  with 
grain- scale  model  resolution  over  regions  experiencing  dynamic  failure. 

1.2  Proposed  Approach 

As  a  precursor  to  a  three-dimensional  (3D)  finite  strain  micromorphic  plasticity  model  (9)  and 
finite  element  (FE)  implementation  (10),  and  overlap  coupling  with  underlying  3D  FE  or  discrete 
element  (DE)  DNS  region,  we  consider  a  simpler,  one  dimensional  (ID)  problem:  overlap 
coupling  between  a  micropolar  linear  elastic  ID  mixed  FE  model  and  a  ID  string  of  Hertzian 
(11,  12)  nonlinear  elastic  DE  spheres. 

To  illustrate  the  application  of  the  micromorphic  plasticity  model  (of  which  micropolar  elasticity 
is  a  subset)  to  the  problem  of  interest,  we  refer  to  an  illustration  in  figure  2  of  a  concurrent 
multiscale  modeling  framework  for  bound  particulate  materials  (target)  impacted  by  a  deformable 
solid  (projectile).  The  higher  order  continuum  micromorphic  plasticity  model  is  used  in  the 
overlap  region  between  a  continuum  finite  element  (FE)  and  DNS  representation  of  the  particulate 
material.  The  additional  degrees  of  freedom  provided  by  the  micromorphic  model  (micro-shear, 
micro-dilation/compaction,  and  micro-rotation)  will  allow  the  overlap  region  to  be  placed  closer 
to  the  region  of  interest,  such  as  at  a  projectile-target  interface.  Further  from  this  interface  region, 
standard  continuum  mechanics  and  constitutive  models  can  be  used.  The  discrete  element  (DE) 
and/or  finite  element  (FE)  representation  of  the  particulate  micro- structure  is  intentionally  not 
shown  in  order  not  to  clutter  the  drawing  of  the  micro-structure.  The  grains  (binder  matrix  not 
shown)  of  the  micro- structure  are  ‘meshed’  using  DEs  and/or  FEs  with  cohesive  surface  elements 
(CSEs).  The  open  circles  denote  continuum  FE  nodes  that  have  prescribed  degrees  of  freedom 
(dofs)  D  based  on  the  underlying  grain-scale  response,  while  the  solid  circles  denote  continuum 
FE  nodes  that  have  free  dofs  D  governed  by  the  micromorphic  continuum  model.  We 

2This  is  a  problem  especially  for  modeling  fragmentation  and  comminution  micro-structurally. 
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intentionally  leave  an  ‘open  window’  (i.e.,  DNS)  on  the  particulate  micro- structural  mesh  in  order 
to  model  dynamic  failure.  If  the  continuum  mesh  overlays  the  whole  particulate  micro- structural 
region,  as  in  (13)  for  atomistic-continuum  coupling,  then  the  continuum  FEs  would  eventually 
become  too  deformed  by  following  the  micro-structural  motion  during  fragmentation.  The 
blue-dashed  box  at  the  bottom-center  of  the  illustration  is  a  micromorphic  continuum  FE  region 
that  can  be  converted  to  a  DNS  region  for  adaptive  high-fidelity  material  modeling  as  the 
projectile  penetrates  the  target. 
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Figure  2.  2D  illustration  of  concurrent  computational  multi-scale  modeling  approach  in  the  contact  inter¬ 
face  region  between  a  bound  particulate  material  (e.g.,  ceramic  target)  and  deformable  solid  body  (e.g., 
refractory  metal  projectile). 


1.3  Focus  of  Report 

Regarding  the  approach  described  in  section  1.2,  this  report  focusses  on  the  ID  overlap  coupling 
between  a  micropolar  linear  elastic  FE  model  and  a  ID  string  of  Hertzian  nonlinear  elastic  DE 
spheres.  An  outline  of  the  report  is  as  follows:  section  2. 1  summarizes  the  Statement  of  Work 
(SOW)  and  the  Tasks,  section  2.2  the  ID  micropolar  linear  elasticity  derived  from  the  3D 
micropolar  theory  using  Timoshenko  beam  kinematics  with  axial  stretch,  section  2.3  the 
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nonlinear  elastic  theory  and  implementation  for  spherical  DE  particles  with  Hertzian  elastic 
contact,  section  2.4  the  overlap  coupling  methodology  and  numerical  examples,  and  section  3. 
summarizes  the  results,  conclusions,  and  future  work. 

The  DE-FE  coupled  implementation  is  currently  limited  to  nonlinear  quasi-statics,  but  the 
formulation  has  been  provided  in  general  for  nonlinear  dynamics,  and  will  be  extended  to 
dynamics  in  future  work. 

1.4  Notation 

Index  notation  will  be  used  wherever  needed  to  clarify  the  presentation.  Cartesian  coordinates 
are  assumed,  so  all  indices  are  subscripts,  and  spatial  partial  derivative  is  the  same  as  covariant 
derivative  (14).  Some  symbolic/direct  notation  is  also  given,  such  that  ( ab)ik  =  atjbjk, 

(a  <g)  b)ijki  =  dijbki.  Boldface  denotes  a  tensor  or  vector,  where  its  index  notation  is  given. 
Subscript  (•)_*  implies  a  spatial  partial  derivative.  Superposed  dot  (□)  =  D(D)/Dt  denotes 

def 

material  time  derivative.  The  symbol  =  implies  a  definition. 
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2.  Technical  Discussion 


2.1  Statement  of  Work  (SOW)  and  Specific  Tasks 

Bound  particulate  materials  are  commonly  found  in  industrial  products,  construction  materials, 
and  nature  (e.g.,  geological  materials).  They  include  poly  crystalline  ceramics  (e.g.,  crystalline 
grains  with  amorphous  grain  boundary  phases),  energetic  materials  (high  explosives  and  solid 
rocket  propellant),  hot  asphalt,  asphalt  pavement  (after  asphalt  has  cured),  mortar,  conventional 
quasi-brittle  concrete,  ductile  fiber  composite  concretes,  and  sandstones,  for  instance.  Bound 
particulate  materials  contain  particles3  (quasi-brittle  or  ductile)  bound  by  binder  material 
oftentimes  called  the  “matrix”. 

The  heterogeneous  nature  of  bound  particulate  materials  governs  its  mechanical  behavior  at  the 
particle-  to  continuum-scales.  The  particle-scale  is  denoted  as  the  scale  at  which  particle-matrix 
mechanical  behavior  is  dominant,  thus  necessitating  that  particles  and  matrix  material  be  resolved 
explicitly  (i.e.,  meshed  directly  in  a  numerical  model),  accounting  for  their  interfaces  and 
differences  in  material  properties.  Currently,  there  is  no  approach  enabling  prediction  of 
initiation  and  propagation  of  dynamic  fracture  in  bound  particulate  materials — for  example, 
poly  crystalline  ceramics,  particulate  energetic  materials,  mortar,  and  sandstone — accounting  for 
their  underlying  particulate  micro  structure  across  multiple  length-scales  concurrently. 

Traditional  continuum  methods  have  provided  the  basis  for  understanding  the  dynamic  fracture  of 
these  materials,  but  cannot  predict  the  initiation  of  dynamic  fracture  without  accounting  for  the 
material’s  particulate  nature.  Direct  numerical  simulation  (DNS)  of  deformation,  intra-particle 
cracking,  and  inter-particle-matrix/binder  debonding  at  the  particle-scale  is  limited  by  current 
computing  power  (even  massively-parallel  computing)  to  a  small  representative  volume  element 
(RVE)  of  the  material,  and  usually  must  assume  overly-restrictive  boundary  conditions  (BCs)  on 
the  RVE  (e.g.,  fixed  normal  displacement). 

Multiscale  modeling  techniques  are  clearly  needed  to  accurately  capture  the  response  of  bound 
particulate  materials  in  a  way  accounting  simultaneously  for  effects  of  the  microstructure  at  the 
particle-scale  and  boundary  conditions  applied  to  the  engineering  structure  of  interest,  at  the 
continuum-scale.  The  services  of  a  scientist  or  engineer  are  required  to  develop  the  mathematical 
theory  and  numerical  methodology  for  multiscale  modeling  of  bound  particulate  materials  of 
interest  to  the  U.S.  Army  Research  Laboratory  (ARL). 

3We  use  ‘particle’  and  ‘grain’  interchangeably. 


6 


The  overall  objective  of  the  proposed  research  is  to  develop  a  concurrent  multi-scale 
computational  modeling  approach  that  couples  regions  of  continuum  deformation  to  regions  of 
particle-matrix  deformation,  cracking,  and  debonding,  while  bridging  the  particle-  to 
continuum-scale  mechanics  to  allow  numerical  adaptivity  in  modeling  initiation  of  dynamic 
fracture  and  degradation  in  bound  particulate  materials. 

For  computational  efficiency,  the  solicited  research  will  use  DNS  only  in  spatial  regions  of 
interest,  such  as  the  initiation  site  of  a  crack  and  its  tip  during  propagation,  and  a  micromorphic 
continuum  approach  will  be  used  in  the  overlap  and  adjacent  regions  to  provide  proper  BCs  on 
the  DNS  region,  as  well  as  an  overlay  continuum  to  which  to  project  the  underlying  particle-scale 
mechanical  response  (stress,  internal  state  variables  (IS Vs)).  The  micromorphic  continuum 
constitutive  model  will  account  for  the  inherent  length  scale  of  damaged  fracture  zone  at  the 
particle-scale,  and,  thus,  includes  the  kinematics  to  enable  the  proper  coupling  with  the  fractured 
DNS  particle  region.  Outside  of  the  DNS  region,  a  micromorphic  extension  of  existing 
continuum  model(s),  with  the  particular  model(s)  to  be  determined  based  on  ARL  needs,  of 
material  behavior  will  be  used. 

This  SOW  calls  for  development  of  the  formulation  and  finite  element  implementation  of  a  finite 
strain  micromorphic  inelastic  constitutive  model  to  bridge  particle-scale  mechanics  to  the 
continuum-scale.  The  desired  result  is  formulation  of  such  a  model,  enabling  a  more  complete 
understanding  of  the  role  of  microstructure-scale  physics  on  the  thermomechanical  properties  and 
performance  of  heterogeneous  materials  of  interest  to  ARL.  These  materials  could  include,  but 
are  not  limited  to,  the  following:  ceramic  materials,  energetic  materials,  geological  materials, 
and  urban  structural  materials. 

2.1.1  Specific  Tasks 

What  follows  is  a  list  of  specific  tasks,  and  a  summary  of  what  was  accomplished  for  each  task. 

1.  Finite  element  implementation  of  finite  strain  micromorphic  pressure-sensitive 
elasto-plasticity  model  (Regueiro,  J.  Eng.  Mech.,  2009)  in  the  continuum  mechanics  code 
Tahoe  source  forge  .  net/projects/tahoe. 

This  implementation  is  ongoing,  but  has  not  been  completed  in  time  for  this  report. 

2.  Interact  with  ARL  researchers  in  order  to  improve  mutual  understanding  with  regards  to 
dynamic  fracture  and  material  degradation  in  heterogeneous  and  particulate  materials  and 
associated  numerical  modeling  techniques. 
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Continue  to  interact  with  ARL  researchers  regarding  their  needs  for  this  research  problem. 

3.  Refine  formulation  of  algorithm  to  couple  finite  strain  micromorphic  continuum  finite 
elements  to  DNS  finite  elements  of  bound  particulate  material  through  an  overlapping 
region. 

This  has  been  accomplished  to  some  extent  through  the  ID  overlap  coupling  problem 
described  in  this  report  in  sections  2.2,  2.3,  and  2.4. 

4.  Implement  coupling  algorithm  in  using  finite  element  code  Tahoe  (both  for  micromorphic 
continuum  and  DNS).  Extension  can  be  made  for  coupling  micromorphic  model  (Tahoe)  to 
DNS  model  (ARL  or  other  finite  element,  or  particle/meshfree,  code). 

The  coupling  algorithm  has  been  implemented  in  a  separate  standalone  code  instead  of 
Tahoe.  Future  research  will  establish  the  coupling  in  Tahoe,  as  part  of  currently  funded 
research  projects. 

2.2  One-dimensional  (ID)  Micropolar  Linear  Isotropic  Elasticity 

This  section  (in  section  2.2.1)  briefly  presents  the  three-dimensional  (3D),  small  strain,  linear 
isotropic  micropolar  elasticity  model  and  balance  of  linear  and  angular  momentum  equations 
based  on  the  work  of  (15),  and  then  provides  more  details  on  the  reduction  to  a  ID  form  (section 
2.2.2)  using  Timoshenko  beam  kinematics  with  axial  stretch  (16).  Finally,  in  section  2.2.3,  the 
ID  form  of  the  model  is  expressed  in  weak  and  Galerkin  forms  for  finite  element  implementation, 
a  mixed  ID  element  is  used  to  interpolate  the  fields,  and  a  numerical  example  is  presented  to 
demonstrate  convergence  of  the  FE  implementation. 

2.2.1  Three-dimensional  (3D)  Micropolar  Linear  Isotropic  Elasticity  and  Balance 


Equations 


The  balance  of  linear  and  angular  momentum,  respectively,  for  a  micropolar  continuum  are  (15) 


&ik,i  +  pbk  ~  puk  —  0 
'mik,l  T  &kmn@ mn  T  pK-k  Pfik  0 


(1) 

(2) 


where  cr/fc  is  the  unsymmetric  Cauchy  stress  tensor  over  body  B,  p  is  the  mass  density,  6/.  is  a  body 
force  per  unit  mass,  Uk  is  the  displacement  vector,  /'//,  is  the  acceleration  vector,  m/fc  is  the 
unsymmetric  couple  stress,  ekmn  is  the  permutation  operator,  £k  is  the  body  couple  per  unit  mass, 
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(5k  =  jfik  is  the  intrinsic  spin  per  unit  mass,  j  is  the  spin  inertia  for  a  spin-isotropic  material  (77), 
ipi  is  the  micro-rotation  vector,  indices  k,  /,-■•  =  1,  2,  3,  and  («))Z  =  d(»)/dxi  denotes  partial 
differentiation  with  respect  to  the  coordinate  xi. 

For  a  linear  isotropic  elastic  micropolar  solid,  the  constitutive  equations  are  (15) 


&kl  Acr?.(5fcZ  T  (2/i  T  T  ^^klmij'm  ^Pm)  (3) 

mki  =  a(pr>rSki  +  (3tpk,i  +  lP>i,k  (4) 


where  A,  /x,  k,  a,  (5,  and  7  are  isotropic  elastic  parameters,  and  the  deformations  are 

ekl  =  2  (Uk,l  +  Ul,k)  (5) 

T k  2  *~’klm'UJm,l  (6) 

Tkl  2^Uk’l  ^l,k)  ® klmT  m  (7) 

where  eki  is  the  classical  small  strain  tensor,  rk  the  axial  vector,  and  rki  the  rotation  tensor. 

2.2.2  ID  Timoshenko  Beam  Kinematics  with  Axial  Stretch  and  Resulting  ID  Micropolar 
Linear  Elasticity 

We  consider  Timoshenko  beam  kinematics  from  (16),  with  superimposed  axial  stretch  for  small 
deformations.  The  displacement  vector  u  and  micro-rotation  vector  ip  are 


Ui 

u  —  x2  0 

<Pi 

'  0  ' 

u  = 

U2 

= 

V 

,  <P  = 

P>2 

= 

0 

.  U3  . 

0 

.  ¥>3  . 

_  9  _ 

where  u  is  the  stretch  in  the  x\  direction,  9  is  the  rotation  of  the  centroidal  axis  about  the  x3  axis, 
v  is  the  transverse  displacement  in  the  x2  direction,  we  ignore  displacement  w  ~  0  in  the  x3 

def 

direction,  and  we  assume  the  micro-rotation  tp3  =  9. 
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Taking  spatial  derivatives,  we  have  the  following  deformation  measures: 


[Uk,l]  — 


p--x2-^~  -e  o 

ox  i  ^  ox  i 


dv 

9xi 

0 


0 

0 


0 

0 


,  {tk,l\  = 


[x k]  0  [&klm'U'm,l\ 


kl\  \&klmXm\ 


\{ik  +  d 


U£k  +  6 


[CfcZ?n^m]  — 


0  0  0 

-0  0  0 

0  0  0 


du  89  ]_  (  dv  n 

dx i  X2dx1  2  l  ° 


1  (  dv  n 

2  9xi  “ 

0 


0 

0 


-Hfe  +  0 


0 

0 

0 


0 

0 


0 

0 

0 


(9) 


(10) 


(11) 


The  unsymmetric  stress  tensor  components  then  result  as 


CTll 

022 

033 

0^12 

021 

mi3 

m3 1 


(A  +  2/i  +  k)cii 

Acn 

Aen 

(2/1  +  K)ei2  -  «Xi2  -  KeUmVm 

(2/1  +  «)e2i  -  nr21  -  Ke21 m</?m 


7' 

0 


90 

90 

9xi 


(/!  +  «)( 


9u 

9xi 


/17 


sh 


f) 


(/i  +  «)7sh 


(12) 

(13) 

(14) 

(15) 

(16) 

(17) 

(18) 


where  7sh  =f  ^  -  0,  cr23  =  0-32  =  0-13  =  ^31  =  0,  and 
mn  =  m22  =  m33  =  m12  =  m21  =  m23  =  m32  =  0. 
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The  balance  of  linear  and  angular  momentum  equations  can  likewise  be  reduced  as 


011,1  -  pui  =  0 
0" 12,1  —  PU  2  =  0 
TO13.1  +  ^12  -  0*21  ~  P$ 3  =  0 


(19) 

(20) 
(21) 


where  $3  =  jp>3  =  jd,  assuming  spin-inertia  j  is  constant  at  small  strains.  However,  to  reach  a 
form  amenable  to  a  ID  mixed  FE  formulation,  we  express  the  weak  form  in  terms  of  the  reduced 
kinematics.  Consider  the  weak  form  in  3D  as 


/  pwkukdv  + 

Jb  JB 

/  PVkftkdv  +  /  rjkjmikdv  —  /  pkek 

Ib  Jb  Jb 


where  the  traction  tk  =  aiknt,  surface  couple  stress  rk 


/  Wktkda 

(22) 

JTt 

mn^mndv  /  TjkVfcdCL 

JTr 

(23) 

=  mikni,  and  the  weighting  functions  are 


W\ 

5u\ 

5u  - 

-  X2  89 

[Wk]  = 

w2 

= 

5u2 

= 

5v 

0 

0 

0 

'  0  ' 

0 

0 

[Vk\  = 

0 

= 

0 

= 

0 

-  93  . 

.  S(P3  . 

.  66  . 

(24) 


(25) 


where  <$(•)  is  used  here  as  a  variational  operator.  Likewise,  u\  =  u  —  xj)  and  u2  =  v.  We 
analyze  each  term  separately,  such  that  the  first  term  of  equation  22  is 


pwkukdv  =  /  p(wiUi  +  w-zu^dv 


' B 


P 


5u(u  —  x2  9)  +  8v(v)  +  59(—X2U  +  (x2  )20) 


dv 


(26) 

(27) 


where  we  consider  that  all  variables  are  functions  only  along  the  1-D  length  xx  (which  we  will 
simplify  as  x). 


11 


Thus,  when  reducing  the  integral  JB(»)dv  =  f L  fA(»)dadx,  we  have 


'  L  L 


pwkUkdv  =  /  5u(pAu  —  pQO)  +  5v(pAv)  +  56(—pQu  +  pIO ) 


dx 


(28) 


where,  the  first  moment  about  the  axis  Q  and  moment  of  inertia  about  Xj  axis  /  are  defined  as 


j  def 


(: x2)2da 


(29) 


Likewise,  the  stress  term  is 


wk,io-ikdv  =  /  (wltlan  +  w2jla12)dv 


(30) 


(A  T  2p  -\-  K^i^Au^x  QO^x)  dv^x^P  Aip) 


+<50ia;(A  +  2/i  +  n)(y—Qu)X  +  /((i)]  dx 


(31) 


Considering  that  the  traction  acts  in  the  x\  direction  on  the  xA  face,  then  the  unit  normal  to  T,  is 

n  —  [1  0  0]r,  and  ti  =  an,  f2  =  0"i2-  Then, 


wktkda  =  /  (witi  +  w2t2)da 


'rt 


Suanda  —  /  58x2anda  +  /  Sva^da 


I A 


(32) 

(33) 


where,  if  we  consider  the  traction  acts  at  x  =  L,  then  concentrated  axial  force  FL  =  JA  anda,  end 
moment  —  ML  =  fA  x2anda,  and  end  shear  force  Vl  =  JA  cr12dA,  and,  thus, 


'UJktkda 


>T  t 


5ulFl  +  56lMl  +  SvLVL 


(34) 


For  the  balance  of  angular  momentum  weak  form,  we  have  for  the  micro-inertia  term 


PVkPkdv  =  /  58(pAjd)dx 


(35) 
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and  for  the  couple  stress  term 


/  V kjmikdv  =  /  r]3tlm13dv  =  /  50.X(~/A0.x)dx 

'b  Jb  Jl 


(36) 


and  the  skew  part  of  the  stress 


-  /  r}kekmn(Jmndv  =  /  86(<Ji2  -  a2i )dv  =  59  [kA{6  -  dx  (37) 

Jb  Jb  Jl 


where  if  the  traction  couple  stress  acts  on  the  X\  face  at  x  =  L,  then 


r)krkda  = 


'rr 


/  59rrii3da  =  59LMeL  ,  =  /  m13<ia 

'  A  J  A 


(38) 


Next,  we  put  these  equations  in  Finite  Element  (FE)  matrix  form. 


2.2.3  Finite  Element  (FE)  Implementation  of  ID  Micropolar  Linear  Elasticity 

The  Galerkin  form  of  the  reduced  equations  may  be  written  as 


5uh(pAuh  -  pQ9h)  +  5vh(pAvh )  +  59h{-pQuh  +  pI6h) 

5uhx{\  +  2/r  +  n)(Auhx  -  Q6hx )  +  5vhx(p  +  K)A{vhx  -  6h) 

+S6hx( A  +  2/r  +  n)(-Quhx  +  I9hx )]  dx  =  5uhLFL  +  59hLML  +  8vhLVL 


'  L  L 

.h 


89h(pAjOh )  +  59hx(^A9hx)  +  [/tA(0h  -  vhx )] 


dx  =  56hLMeL 


(39) 

(40) 
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The  interpolations,  and  their  spatial  derivatives,  for  a  mixed  element  shown  in  figure  3  are  the 
following 


'x(a) 


a=  1 
2 


7V“  N% 


rle 

ax{l) 

r]e 

ax(2) 


=  Nu’edex 


e 

y(a ) 


N X  N%  N% 


e\i)  = 


Suh(0 

Svh(0 

69h(0 


a=  1 
2 

a=  1 

Nu’eci 


Nv,eCy 


ay(  i) 

rle 

ay(  2) 
r\e 

ay{  3) 


=  Nv,edey 


6»(a) 


N°  N° 


de 

u0(i) 

r\e 

ae(2)  J 


=  Nd’edeg 


=  Nd’ec 


«*(£),» 

=  Buedex 

Bue  = 

1 - 

AO,* 

=  Bv,edy  , 

Bve  = 

_Nvhx  Nlx 

=  B^dl, 

Be,e  = 

ni 

Suh(0,x 

=  Bu,ecex 

Sv\0,x 

=  Bv,eCy 

seh(0,x 

=  B9’eceg 

K 


3,x 


ef 


l® 


vhe 
n  v3 


he 


2  '  0? 


Figure  3.  Finite  element  degrees  of  freedom  (dof)  for  mixed  formulation  Timoshenko  beam  with  axial 
stretch.  The  middle  node  3  is  at  the  center  of  the  element,  i.e.,  at  £  =  0,  where  £  is  the  natural  coordinate 
(3). 


14 


We  assign  N0,e  =  Nu,e  and  B9,e  =  Bu,e,  and  the  element  dofs  are 


1 

"Ci  rH 

1 - 

< 

de 

Ux(2) 

< 

de 

ay{  i) 

h  e 

V 2 

=  de  = 

de 

ay(  2) 

he 

de 

ay(  3) 

Of 

de 

u0(i) 

.  . 

de 

L  a0(2)  J 

dex 

dl 


(41) 


Substituting  these  interpolations  into  the  Galerkin  form  of  balance  of  linear  momentum  in 
equation  39,  and  grouping  terms,  we  can  define  element  mass  and  stiffness  matrices  as 


ne\ 


e=l 


pAe(Nu,e)T  Nu,edx^j  dex  -  pQe(Nu,e)T  Nd,edx^j 


6 


V 

ryy-yiiu,e. 


mu 


+  (jf  (A  +  2p  +  K)Ae(Bu’e)T Bu’edx^J  dex  -  (A  +  2//  +  n)Qe{Bu’e)T Be'ed^j  de9 


k? 


k 


u9,e 


=  /j 


(42) 


"el 


e\T 


K) 


e=l 


pAe(Nv,e)T Nv,edx^j  dey  +  Qf  (/i  +  K)Ae(Bv’e)TBv’edx^J  d 


v 


k 


vv.e 


(/i  +  n)Ae(Bv'e)T Ne’edx  )ded  =  f 


k 


v9,e 


(43) 


™el 


A(c*) 


,e\T 


e=l 


pQe(Ne,e)T Nu,edx^j  dex  +  p/e(AT4e)TAr0’e(ix^  d; 


••,e 

0 


m" 


m" 


-  (jf  (A  +  2/r  +  K)Qe(B0’e)TBu'e<hJ  +  (jf  (A  +  2/r  +  K,)Ie(Be’e)T B0>edfoJ  dee 


j^6u,e 


k°ei,e 


=  f 


Ml 


(44) 
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where  A 's  the  element  assembly  operator,  ne\  the  number  of  elements,  and  te  =  he  the  length  of 
an  element  e,  where  the  equations  can  be  written  more  concisely  in  matrix  form  as 


ree  i 


Alc 

e=l 
riel 


S)T 

e\T 


muu,ecTx  -  mud’ede9  +  kuu’edex  -  kuU'ede9  =  f 


uQ.e  je 


e=l 

riel 


e=l 


j\(cey)T  mvv,eidy  +  kvv’edey  -  kv9’edee  =  fv 


„e\T 


-m9u'edex  +  m9ei’eded  -  k9u’edex  +  km’ede6  =  f 


(45) 

(46) 

(47) 


Likewise,  when  substituting  the  interpolations  into  the  Galerkin  form  of  balance  of  angular 
momentum  in  equation  40,  and  grouping  terms,  we  have 


«el 


e=l 


~e\T 


pje(N9,e)T Ne,edx^j  de9  +  ^Ae{B9’e)TB9’edx^j  dee 


rrv< 


k 


992, t 


+  KAe(N9’e)TN9’edx ^  deg  -  Qf  K,Ae(N9’e)T Bv,edx'SJ  dey  =  fMe 


k 


993,  e 


fc9 


(48) 


where,  in  summary,  we  have 


»»el 


e=l 


+  fcW2’e4  +  k993’edee  -  k9v'edey  =  f 


Me 


(49) 


Adding  equations  47  and  49  for  the  micro-rotation  dofs,  we  have 


reel 


e—1 


~e\T 


-m9u’edex  +  ( m991'e  +  m992’e)deg  -  k9u’edex  +  ( k991’e  +  k992'e  +  k993’e)ded 
~k9v’edey  =  feM  +  />]  (50) 
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Accounting  for  essential  boundary  conditions  (BCs),  and  assembling  the  global  FE  matrix 
equations  (3),  we  arrive  at  the  coupled  system  of  matrix  FE  equations  to  solve  for  the  unknown 
dofs  as 


MdD  +  KdD  =  Fd  (51) 


dx 

FF 

D  = 

dy 

,  D  = 

dy 

,  Fd  = 

Fv 

dg 

dg 

Fm  +  F  Me 

Muu 

0 

-Mue 

0 

Mvv 

0 

- M9u 

0 

M" 

J£ILU 

0 

-Kue 

0 

J£VV 

—Kv0 

-k9u 

-K0v 

K" 

Since  Nd’e  =  Nu'e,  then  M9u  =  Mu 6  and  MD  is  symmetric,  and  since  Be’e  =  Bu  e  then 
K9u  =  Ku9 ,  but  Kd  is,  in  general,  unsymmetric  because  fi  >  0  which  leads  to  K9v  ^  KL'e . 

Given  the  mixed  ID  Timoshenko  beam  micropolar  elastic  finite  element  with  axial  stretch  in 
figure  3,  we  can  select  specific  shape  functions,  with  resulting  first  spatial  derivatives,  as  follows: 


Nu’e  = 

Nv,e  = 


1  r 

2 


TDIL  ,e 

B  =  ~he 


1-e  1  +  ^ 

K(e-l)  |^+1)  l~e2 


-1  1 


(52) 

(53) 


where  he  is  the  element  length. 
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With  these  element  interpolation  matrices  Ne  and  ‘strain-displacement’  matrices  B'  ,  and  a  bar 
with  circular  cross-section  such  that  the  first  moment  Qe  =  0  (see  section  2.2.4),  we  arrive  at 
specific  forms  of  the  element  stiffness  matrices  as 


fouu,e  _ 

kvv,e  = 

kv6,e  = 

k881,e  = 
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—(A  +  2/i  +  k) 
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_  5 
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0 
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6  12 


(54) 

(55) 

(56) 

(57) 

(58) 


We  can  similarly  derive  the  element  mass  matrices,  but  since  our  simulations  are  currently  limited 
to  quasi-static  problems,  we  do  not  show  all  details  for  including  the  inertia  terms.  This  is  part  of 
future  work. 

2.2.4  Convergence  of  ID  Micropolar  Linear  Elastic  FE 

In  this  section,  we  take  the  FE  formulation  and  implementation  from  the  previous  section  and  test 
its  convergence  with  regard  to  spatial  discretization  refinement  (i.e.,  he  — >  0).  It  is  well  known 
that  in  the  thin  limit  the  Timoshenko  beam  formulation  will  do  a  poor  job  calculating  transverse 
displacement  (16),  because  the  classical  equal  interpolation  (Nu,e  =  Nv,e,  and  Bu,e  =  Bv,e  ) 
element  “locks”  as  the  transverse  shear  strain  approaches  zero  in  the  thin  limit.  Thus,  we  use  the 
mixed  formulation  (Nu,e  Nv,e,  and  B"  f  ^  Bv  e )  in  figure  3  and  reduced  integration  to 
alleviate  this  problem4  (16).  Figure  4  shows  the  five-element  mesh,  and  the  force  versus 
displacement  results  for  two-,  five-,  ten-,  and  twenty-element  meshes.  We  can  see  that  for  the 
axial  force  and  displacement,  the  result  is  the  same  for  all  meshes  (as  expected  for  small  strain 
theory,  where  axial  and  transverse  displacements  are  decoupled),  whereas  for  the  five-element 

4We  note  that  in  our  case  the  reduced  integration  gives  us  no  added  benefit,  as  the  mixed  formulation  seems 
sufficient  to  address  any  potential  locking.  See  results  in  figure  4. 
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mesh,  for  the  transverse  force  and  displacement,  the  results  appear  convergent.  In  the  overlap 
coupling  simulations  in  section  2.4,  we  will  use  the  five-element  mesh. 


Each  mesh  is  L  =  20cm  in  length,  with  rod  circular  cross-section  with  radius  R  =  0.5cm,  as 
shown  in  figure  4.  The  applied  forces  are  F^X1  =  lkN,  FyXI  =  lkN.  For  the  elastic 
parameters,  we  use  Young’s  modulus  E  =  lOGPa,  Poisson’s  ratio  v  =  0.3  (approximate  for 
quartz),  and  approximate  A  and  /i  as 


A 


Ev 

(l  +  z/)(l-2 v)  ’ 


E 

2(l  +  i/) 


(59) 


even  though  A  and  //  are  not  the  Lame  parameters.  We  approximate  k  =  0.1  fi,  and 
7  =  (0.05 )P/i,  with  elastic  length  scale  i  =  1cm.  The  cross-sectional  area  A  =  nR2  =  7.85e-5 
m2,  first  moment  about  the  x-axis  Q  —  0,  and  moment  of  inertia  about  the  x-axis 
I  =  jttR4  =  4.9e-10  m4. 


pEXT 
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f  vs  d  of  node  1 

X  X 


20cm  - 

- - ©■ 


■9- 
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f  vs  d  of  node  1 

y  y 


Figure  4.  (top)  Five-eiement  mesh,  (bottom)  Demonstration  of  convergence  of  thin  Timoshenko  beam 
mixed  FE  implementation  in  compression,  bending,  and  shear.  Axial  force  versus  displacement  of  the  left 
end  node  is  exact,  whereas  the  two-element  mesh  may  be  too  coarse,  while  the  live-,  ten-,  and  twenty- 
element  meshes  give  nearly  the  same  transverse  displacement  upon  an  applied  transverse  force  FyXT. 
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2.3  ID  String  of  Hertzian  Nonlinear  Elastic  Discrete  Element  (DE)  Spheres 

We  refer  to  figure  5  for  the  kinematics  and  forces  at  contact  e  between  two  discrete  element  (DE) 
spheres  a  and  f3.  From  Hertz-Mindlin  elastic  contact  theory  between  stiff  spheres  (11,  12,  18), 
we  have 
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i£0^/2K)3/2 

E°  =  9(1-^)  5  R°  =  R /2  ’  Sn  =  (lx  ~  € 
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y  =  57^7 ,  «  =  K)1/2 ,  S',  =  <  -  qi  +  R(oia  +  (J) 
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(61) 


where  E  is  the  Young’s  modulus,  u  is  the  Poisson’s  ratio,  and  R  is  the  radius  of  a  spherical 
particle  (particles  are  assumed  to  have  equal  R  in  this  case).  We  can  then  assemble  the  internal 
force  vector  fINT’e  and  local  consistent  tangent  dfINT,e/dqe  associated  with  dof  vector  qf  at 
contact  e  as  follows: 
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where  the  contact  moment  ff  is  calculated  by  factoring  a  rotational  stiffness  with  dimensionless 
scalar  B,  such  that 

,fl  =  B^R5/2(coa  +  c /)  +  Rfy  (66) 

This  is  done  to  avoid  a  rank  deficient  local  consistent  tangent  df1NT’e/dqe  because  of  the  linear 
dependence  of  Rfy  on  fy. 


Figure  5.  Kinematics  and  forces  of  two  DE  spheres  a  and  (3  contacting  at  contact  e. 

These  internal  force  and  moment  vectors  are  assembled  into  a  global  nonlinear  internal  force  and 
moment  vector  F  ( q ),  that  when  combined  with  an  external  force  and  moment  vector  F  , 

lead  to  a  residual  form  of  the  balance  of  linear  and  angular  momentum  to  solve  using  the 
Newton-Raphson  method, 
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(67) 


For  the  external  force  and  moment  vector  Fl'  x  r ,  we  will  insert  the  boundary  conditions  directly 
into  the  corresponding  global  dofs  in  FhX1 .  Simulations  using  the  DE  implementation  will  be 
demonstrated  in  the  context  of  the  overlap  coupling  methodology  discussed  in  the  next  section. 
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2.4  Overlap  Coupling  Between  ID  Micropolar  FEs  and  a  String  of  Spherical  DEs 

An  aspect  of  the  computational  concurrent  multiscale  modeling  approach  is  to  couple  regions  of 
material  represented  by  particles,  Discrete  Element  (DE),  to  regions  of  material  represented  by 
continuum.  Finite  Element  (EE).  Another  aspect  is  to  bridge  the  particle  mechanics  to  a 
continuum  representation  using  finite  strain  micromorphic  elasto-plasticity,  whereas  the  small 
strain  micropolar  continuum  is  a  simple  approximation  of  stiff  particles  with  small  frictional 
sliding  in  the  overlap  region  (we  consider  no  sliding  in  the  numerical  examples  in  sections  2.4.3 
and  2.4.4).  The  coupling  implementation  will  allow  arbitrarily  overlapping  particle  and 
continuum  regions  in  a  single  “hand- shaking”  or  overlap  region,  such  that  fictitious  forces  and 
wave  reflections  are  minimized  in  the  overlap  region.  In  theory,  for  nearly  homogeneous 
deformation,  if  the  particle  and  continuum  regions  share  the  same  region  (i.e.,  are  completely 
overlapped),  the  results  should  be  the  same  as  if  the  overlap  region  is  a  subset  of  the  overall 
problem  domain  (cf.  figure  6).  This  will  serve  as  a  benchmark  problem  for  the  numerical 
implementation.  The  coupling  implementation  extends  to  particle  mechanics  and  micropolar 
continuum  the  “bridging  scale  decomposition”  proposed  by  (19)  and  modifications  thereof  by 
(13)  (see  references  therein  for  further  background  on  these  methods). 

2.4.1  3D  Kinematics 

Here,  a  summary  of  the  kinematics  of  the  coupled  regions  is  given  for  general  3-D  kinematics, 
following  the  illustration  shown  in  figure  6.  It  is  assumed  that  the  finite  element  mesh  covers  the 
domain  of  the  problem  in  which  the  material  is  behaving  more  solid-like,  whereas  in  regions  of 
large  relative  particle  motion  (fluid-like),  a  particle  mechanics  representation  can  be  used  (in  this 
case,  DE).  In  figure  6,  discrete  domains  are  defined,  where  the  purple  background  denotes  the  FE 
overlap  region  Bh  with  underlying  ghost  particles,  aqua  blue  the  FE  continuum  region  Bh  with  no 
underlying  particles,  and  white  background  (with  brown  particles)  the  free  particle  region 
Bh  U  BDE .  In  summary,  the  finite  element  domain  Bh  is  the  union  of  pure  continuum  FE  domain 
Bh,  overlapping  FE  domain  with  underlying  ghost  particles  Bh,  and  overlapping  FE  domain  with 
underlying  free  particles  Bh,  such  that  Bh  =  Bh  U  Bh  U  Bh .  The  pure  particle  domain  with  no 
overlapping  FE  domain  (i.e.,  the  ‘open-window’)  is  indicated  by  BDE .  The  goal  is  to  have  the 
overlap  region  Bh  U  Bh  as  close  to  the  region  of  interest  (e.g.,  penetrator  skin)  as  to  minimize  the 
number  of  particles,  and,  thus,  computational  effort.  Following  some  of  the  same  notation 
presented  in  (13),  we  define  a  generalized  dof  vector  Q  for  particle  displacements  and  rotations  in 
the  system  as 


Q  =  [ Qa ,  Qp,  ■  ■  ■ ,  <?7,  Up, . . . ,  w7]r,  a,p,...,'yeA 


(68) 
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Q  o  free  particles 
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by  continuum  displacement  and  rotation  fields) 


D%  finite  element  nodes  whose  motion  is  unprescribed 

jj  O  finite  element  nodes  whose  motion  is  prescribed 
by  underlying  particles 


Figure  6.  Two-dimensional  illustration  of  the  coupling  between  particle  and  continuum  regions. 

where  qa  is  the  displacement  vector  of  particle  a,  u)a  its  rotation  vector,  and  A  is  the  set  of  all 
particles.  Likewise,  the  finite  element  nodal  displacements  and  rotations  are  written  as 


D  [dai  ■  ■  ■  j  dCJ  0 0ei  ■  ■  ■  i  0 f\ 

a,b,  ■  ■  ■  ,c  £  Kf  ,  d,  e, . . . ,  f  e  M 


(69) 


where  da  is  the  displacement  vector  of  node  a,  0,i  is  the  rotation  vector  of  node  d,  J\T  is  the  set  of 
all  nodes,  and  A4  is  the  set  of  finite  element  nodes  with  rotational  degrees  of  freedom,  where 
A4  c  J\T.  In  order  to  satisfy  the  boundary  conditions  for  both  regions,  the  motion  of  the  particles 
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in  the  overlap  region  (referred  to  as  “ghost  particles,”  cf.  figure  6)  is  prescribed  by  the  continuum 
displacement  and  rotation  fields,  and  written  as 


Q  =  [qa,q, 3, 


/?)•••)  w. y]T ,  a,  p, . . . ,  7  e  A ,  AeBh 


(70) 


while  the  unprescribed  (or  free)  particle  displacements  and  rotations  are 

Q  =  ■■■,vr]]T,  5,e, ...,qeA,  AeJ3huBDE  (71) 

where  *4.  U  A  =  A  and  A  D  A  =  0.  Likewise,  the  displacements  and  rotations  of  nodes 
overlaying  the  particle  region  are  prescribed  by  the  particle  motion  and  written  as 


D  =  [da,db,...,dc,0d,de,...,0f]T  (72) 

a,  b, . . . ,  c  e  Af  ,  d,  e, . . . ,  /  G  M 

St,  M  e  Bh  u  Bh 


while  the  unprescribed  (or  free)  nodal  displacements  and  rotations  are 


D  =  [dm,dn,...,ds,9t,0u,...,0v]T  (73) 

m,  n, . . . ,  s  G  Af  ,  t ,  u, . . . ,  v  G  A4 
Af,  M  e  Bh  u  Bh 

where  AfL)Af  =  Af,Afr\Af  =  ®,  A4UA4  =  A4,  and  M  fl  M  =  0.  Referring  to  figure  6,  the 
prescribed  particle  motions  Q  can  be  viewed  as  boundary  constraints  on  the  free  particle  region, 
and  likewise  the  prescribed  finite  element  nodal  displacements  and  rotations  D  can  be  viewed  as 
boundary  constraints  on  the  finite  element  mesh  in  the  overlap  region. 

In  general,  the  displacement  vector  of  a  particle  a  can  be  represented  by  the  finite  element 
interpolation  of  the  continuum  macro-displacement  field  uh  evaluated  at  the  particle  centroid  xa, 
such  that 


uh(xa,t )  =  ^2N^(xa)da(t)  ,  a  e  A  (74) 

aeAf 
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where  N are  the  shape  functions  associated  with  the  continuum  displacement  field  uh.  Recall 
that  iV“  have  compact  support  and,  thus,  are  only  evaluated  for  particles  with  centroids  that  lie 
within  an  element  containing  node  a  in  its  domain.  In  DE,  particle  dofs  (translations  and 
rotations)  are  tracked  at  the  particle  centroids,  as  are  resultant  forces  and  moments  (from  forces 
acting  at  contacts).  For  example,  we  can  write  the  prescribed  displacement  of  ghost  particle  a  as 

qa(t)  =  uh(xa,t)  =  22  N%(xa)da(t)  ,  a  e  A  (75) 

o-GA/" 

Likewise,  particle  rotation  vectors  can  be  represented  by  the  finite  element  interpolation  of  the 
continuum  micro-rotation  field  iph  evaluated  at  the  particle  centroid  xa,  such  that 

tph(xQ1t)  =  22  Nb(xa)Ob{t)  ,  aeA  (76) 

beM 

where  N£  are  the  shape  functions  associated  with  the  micro-rotation  field  cph.  For  example,  we 
can  write  the  prescribed  rotation  of  ghost  particle  a  as 

uja{t)  =  Iph(xa,t )  =  22  Nb(.Xa)Ob(t)  ,  a  e  A  (77) 

beM 

For  all  ghost  particles  (cf.  figure  6),  the  interpolations  can  be  written  as 

Q  =  NqD  ■  D  +  NqQ  ■  D  (78) 

where  NqD  and  N qq  are  shape  function  matrices  containing  individual  nodal  shape  functions 
iV“  and  ,  but  for  now  these  matrices  will  be  kept  general  to  increase  our  flexibility  in  choosing 
interpolation/projection  functions  (such  as  those  used  in  meshfree  methods).  Overall,  the  particle 
displacements  and  rotations  may  be  written  as 
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where  Q'  is  introduced  (13)  as  the  error  (or  “fine-scale”  (19))  in  the  interpolation  of  the  free 
particle  displacements  and  rotations  Q,  whose  function  space  is  not  rich  enough  to  represent  the 
true  free  particle  motion.  The  shape  function  matrices  N  are,  in  general,  not  square  because  the 
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number  of  free  particles  are  not  the  same  as  free  nodes  and  prescribed  nodes,  and  the  number  of 
ghost  particles  is  not  the  same  as  prescribed  and  free  nodes.  A  scalar  measure  of  error  in  particle 
displacements  and  rotations  is  defined  as  (13) 


e  =  Q'  Q'  (80) 

which  may  be  minimized  with  respect  to  prescribed  continuum  nodal  dofs  D  to  solve  for  D  in 
terms  of  free  particle  and  continuum  nodal  dofs  as 


D  =  MBBNId(Q  -  N^D)  .  Moo  =  NtqBNqB 


(81) 


This  is  known  as  the  “discretized  L2  projection”  (13)  of  the  free  particle  motion  Q  and  free  nodal 
dofs  D  onto  the  prescribed  nodals  dofs  D.  Upon  substituting  equation  81  into  equation  78,  we 
may  write  the  prescribed  particle  dofs  Q  in  terms  of  free  particle  Q  and  continuum  nodal  D  dofs. 
In  summary,  these  relations  are  written  as 


Q  —  B^qQ  +  B^dD  (82) 

D  =  BfjgQ  +  B5dD  (83) 

where 


Bqq 

—  NqqBQq 

(84) 

BQD 

—  NqD  +  N$5BSd 

(85) 

Bdq 

= 

DD  QD 

(86) 

B5d 

=  -M-±Nt„~Nod 

DD  QD  ^ 

(87) 

As  shown  in  figure  6,  for  a  finite  element  implementation  of  this  dof  coupling,  we  expect  that  free 
particle  dofs  Q  will  not  fall  within  the  support  of  free  continuum  nodal  dofs  D,  such  that  it  can  be 
assumed  that  Nqd  =  0  and 


Q  —  BqqQ  +  BqDD  ,  D  —  B5qQ 


(88) 
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where 


(89) 

(90) 


The  assumption  NqD  ^  0  would  be  valid  for  a  meshfree  projection  of  the  particle  motions  to  the 
FE  nodal  dofs,  as  in  (13),  where  we  could  imagine  that  the  domain  of  influence  of  the  meshfree 
projection  could  encompass  a  free  particle  centroid;  the  degree  of  encompassment  would  be 
controlled  by  the  chosen  support  size  of  the  meshfree  kernel  function.  The  choice  of  meshfree 
projection  in  (13)  was  not  necessarily  to  allow  Q  be  projected  to  D  (and  vice  versa),  but  to 
remove  the  computationally  costly  calculation  of  the  inverse  AT  ~~  in  equations  82  and  83. 

2.4.2  3D  Kinetic  and  Potential  Energy  Partitioning  and  Coupling 

For  the  particle  DE  equations,  the  kinetic  energy  is  T(- ,  dissipation  function  F®,  and  potential 
energy  UQ ,  such  that 


Tq  =  i QMqQ 
Fq  =  aTQ 


(91) 


where  M®  is  a  mass  and  rotary  inertia  matrix  (not  given  in  section  2.3),  and  fint’Q  is  provided 


from  equation  67.  The  dissipation  function  FQ  is  written  as  a  linear  function  of  the  kinetic 
energy  T(-  with  proportionality  coefficient  a,  which  falls  within  the  class  of  damping  called 
“Rayleigh  damping”  (pg.130  of  (20)).  For  the  micropolar  continuum  FE  equations,  TD  is  the 
kinetic  energy,  FD  the  dissipation  function,  and  UD  the  potential  energy,  such  that 


Fd  =  0 


(92) 
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where  MD  and  FINT,D(S)  =  KD S  come  from  equation  51.  We  assume  the  total  kinetic  and 
potential  energy  and  dissipation  of  the  coupled  particle-continuum  system  may  be  written  as  the 
sum  of  the  energies 


t(q,d) 

=  TQ(Q.Q(Q.i)))  +  TD(b,b(Q)) 

(93) 

U{Q,  D) 

=  Uq(Q.Q(Q,D))+Ud(D.D(Q)) 

(94) 

F{Q,D) 

=  F»(Q,Q(Q,D)) 

(95) 

where  we  have  indicated  the  functional  dependence  of  the  prescribed  particle  motion  and  nodal 
dofs  solely  upon  the  free  particle  motion  and  nodal  dofs  Q  and  D.  respectively.  Note  that  the 
dissipation  function  F  =  F(-  only  applies  for  the  particle  system,  and  only  for  static  problems 
(dynamic  relaxation  DE  simulation).  For  purely  dynamical  problems,  FQ  =  0,  and  there  is  only 
dissipation  in  the  particle  system  if  particles  are  allowed  to  slide  frictionally,  and  the  continuum 
has  plasticity  or  other  inelastic  constitutive  response.  Lagrange’s  equations  may  then  be  stated  as 


d  fdT\  _  dT  dF  dU 

dt  V  Dq)  ~  dQ  +  dQ  +  dQ 
d  (dT\  dT  dF  dU 
~dt  \dD )  dD  +  ~dD  +  dD 
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(96) 


which  lead  to  a  coupled  system  of  governing  equations  (linear  and  angular  momentum)  for  the 
coupled  particle-continuum  mechanics.  The  derivatives  are 
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If  the  potential  energy  U  is  nonlinear  with  regard  to  particle  frictional  sliding  and  micropolar  (or 
micromorphic)  plasticity,  then  (96)  may  be  integrated  in  time  and  linearized  for  solution  by  the 
Newton-Raphson  method.  The  benefit  of  this  multiscale  method,  as  pointed  out  by  (79),  is  that 
time  steps  are  different  for  the  DE  and  FE  solutions.  A  multiscale  time  stepping  scheme  will 
follow  an  approach  similar  to  (79).  To  complete  (96)  and  identify  an  approach  to  energy 
partitioning,  the  individual  derivatives  may  be  written  as 


d  TQ 
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(103) 

(104) 
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where  superscript  Q  denotes  free  particle  motion  and  0  ghost  particle  motion,  whereas 
superscript  79  denotes  free  nodal  dofs  and  79  prescribed  nodal  dofs.  The  energy  partitioning  will 
be  introduced  through  the  definition  of  these  terms  below.  First,  substitute  equations  103-107 
into  equation  96  to  arrive  at  the  coupled  nonlinear  equations  in  terms  of  Q  and  D  as 

(m«  +  B|0M<5sgQ  +  BtBqMBBBq)  Q 
+  (BIqmSbqd  +  bIqm~Dbdd )  D 

+  (c°  +  BlQ c5bqq)  Q  +  bTqQcSbqDb 
+f'nt^(Q)  +  bIqf^[bSqQ  +  bSdd] 

+Bt5qFint’B[B3qQ  +  BSdD }  =  Fext’q  +  BlQFEXT$  (108) 

{bSdmSbqd  +  bBdmBbdq)  q 

+  ( Md  +  +  bIdmbbBd )  D 

+bTqdF’NTS\BqqQ  +  Bqdb ]  +  BtBdF,nt-B{BBqQ  +  BBdD] 

+fint,d =  fext,d  +  Bt  FEXT,D  (109) 
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where  an  expression  in  brackets  [•]  denotes  the  functional  dependence  of  the  nonlinear  internal 
force  and  moment  vector.  Note  the  projections  through  the  B  matrices  of  the  corresponding 
mass,  rotary  inertia,  and  damping  matrices,  and  forcing  and  moment  vectors.  First,  starting  with 
the  mass  and  rotary  inertia  matrices  for  the  particles,  we  partition  the  kinetic  energy5  as  follows: 


mq 

—  (1  -  f)  m ^  ,  a  G  A  ,  xa  e  Bh 

OL 

(110) 

mq 

=  Mde’q  +  MQ 

(111) 

mde,q  =  A  >  P  e  A  >  Xp  e  BDE 

0 


MQ  =  (1  -  f)  A  mQp  ,  P  e  A  ,  Xp  e  Bh 

0 

where  is  the  mass  and  rotary  inertia  matrix  of  ghost  particles  in  Bh,  MDE,Q  the  mass  and 
rotary  inertia  matrix  of  free  particles  in  BDE,  M  the  mass  and  rotary  inertia  matrix  of  free 
particles  in  Bh,  and  f  is  a  weighting  factor  for  the  kinetic  energy  in  the  overlap  region  Bh  U  Bh . 
For  no  homogenized  continuum  contribution  to  the  kinetic  energy  in  the  overlap  region,  f  =  0, 
and  for  full  continuum  homogenization  of  the  underlying  particle  kinetic  energy,  r  —  1.  In  our 
case,  we  will  consider  the  range  0  <  f  <  1.  Given  that  the  proposed  multiscale  modeling 
framework  is  to  be  used  in  an  adaptive  fashion  in  the  future,  having  an  overlaying  continuum 
homogenization  of  the  particle  response  is  attractive  when  particle  is  converted  to  continuum 
representation,  and  vice  versa  (in  a  statistical  manner6).  For  the  mass  and  rotary  inertia  matrices 
associated  with  the  micropolar  continuum,  we  partition  the  kinetic  energy  as  follows: 

5  and  dissipation  function  through  the  mass-proportional  damping  for  dynamic  relaxation  of  a  static  DE  analysis 

Statistical,  in  the  sense  of  generating  a  particle  representation  from  a  continuum  one,  where  the  underlying  par¬ 
ticle  system  does  not  exist;  converting  from  particle  to  continuum  representation  is  straightforward  given  the  built-in 
homogenization  that  the  micropolar  continuum  possesses 
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(112) 
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(■ mu’e ) 
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M  —  fa  (f  (rnB'e)  +  remB'e) 
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— -D 

M  =  f 


M  0 


0  M 
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- — D  -  D 

M  +M 


M°  = 


Mu  0 
0  Mv 


,  M  =  (m“’e)  ,  M  =  /\  (mip,e) 

e&Bh  e£Bh 


,  M“  =  A  mu,e  ,  M*  =  fa  m^£ 

e&Bh  e£Bh 


(113) 


where  Af  is  the  continuum  mass  and  rotary  inertia  matrix  associated  with  prescribed  nodal  dofs 

~  —  D  ~ 

in  B  ,  M  the  continuum  mass  and  rotary  inertia  matrix  associated  with  free  nodal  dofs  in  B  , 

- D  - D  -—'D 

where  M  and  M  arc  extracted  from  the  total  mass  and  rotary  inertia  matrix  M  in  B1,  with 
superscript  ( • ) ,J  denoting  the  full  mass  and  rotary  inertia  matrix  associated  with  elements  in  Bh , 
(•)  is  a  homogenization  operator,  re  is  the  partitioning  coefficient  of  continuum  kinetic  energy 
associated  with  element  Be  C  Bh.  A  simple  choice  is  a  volume  fraction  re  =  Be,D /Be,  where 
Be,D  =  Be  —  Be,Q;  Be,D  is  the  non-overlapping  continuum  part  of  element  volume  Be  C  Bh,  and 
Bc  °-  is  the  overlapped  ghost  particle  volume  in  the  element  (cf.  figure  6).  For  kinetic  energy 
partitioning,  a  volume  fraction  that  directly  relates  to  mass  and  rotary  inertia  partitioning  seems 

- D 

an  appropriate  choice.  M  is  the  homogenized  continuum  mass  and  rotary  inertia  matrix 
associated  with  prescribed  nodal  dofs  in  Bh;  where  if  f  =  0,  there  is  no  continuum 
homogenization  in  Bh  (i.e.,  all  kinetic  energy  is  due  to  underlying  particles).  MD  is  the 
continuum  mass  and  rotary  inertia  matrix  associated  with  free  nodal  dofs  in  the  pure  continuum 
FE  domain  Bh. 


31 


For  the  potential  energy  (internal  force  and  moment)  and  external  force  partitioning  in  the  particle 
system,  we  write 


pINT,Q 


p*INT,Q 


pEXT,Q 

pEXT,Q 


(1  -  q)  A  flNT’Q  ,*e£Bh 
6 

p\INT,DE,Q  pINT’Q 

FINT,DE,Q  =  A  fINT,Q  Xs  e  BDE 

FINT’Q  =  (1  -  g)  A  fsNT’Q  ,  *s  eBh 
s 

(1  -  q)  A  f?XT,Q  ,xeeBh 

e 

jpEXT,DE,Q  p,EXT’Q 

FEXT,DE,Q  =  ^  fEXT,Q  Xs  e  BDE 

FEXT,Q  =  (1  -  q)  A  fsXT’Q  i*s£B‘ 
s 


(114) 

(115) 


(116) 

(117) 


where  F,NT'°-  is  the  internal  force  and  moment  vector  associated  with  ghost  particle  contacts  in 
Bh  ,  f1nt’de’Q  is  the  internal  force  and  moment  vector  associated  with  free  particle  contacts  in 

^ INT,Q 

B  .  F  is  the  internal  force  and  moment  vector  associated  with  free  particle  contacts  in  B  , 
pEXT.Q  js  1|1C  externai  force  and  moment  vector  associated  with  ghost  particle  contacts  in  Bh  , 
fext,de,q  js  the  external  force  and  moment  vector  associated  with  free  particle  contacts  in  BDE, 

EXT,Q 

F  is  the  external  force  and  moment  vector  associated  with  free  particle  contacts  in  B  ,  and 
q  is  a  weighting  factor  for  the  potential  energy  in  the  overlap  region  Bh  U  Bh.  For  no 
homogenized  continuum  contribution  to  the  potential  energy  in  the  overlap  region,  q  —  0,  and  for 
full  continuum  homogenization  of  the  underlying  particle  potential  energy,  q  =  1.  In  our  case, 
we  will  consider  the  range  0  <  q  <  1.  Note  that  in  (13),  they  chose  q  =  0.  Their  Cauchy-Born 
elastic  constitutive  model  acts  like  a  homogenization  operator  on  the  underlying  atomistic 
response,  but  instead  of  keeping  an  overlain  Cauchy-Born  representation,  the  potential  energy  is 
completely  represented  by  the  underlying  atomistic  response,  except  in  the  overlap  region  Bh 
where  partitioning  occurs.  Note  that  xe  and  x$  denote  positions  of  particle  contacts  for 
calculating  internal  force  and  moment  vectors  in  equations  114-118,  whereas  xa  and  Xs  in 
equations  1 10  and  111  denote  particle  centroids  for  calculating  particle  mass  and  rotary  inertia 
matrices. 
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For  the  potential  energy  (internal  force)  partitioning  in  the  continuum,  we  write 


pINT,D 


~INT,D  pINT,D 


jINT,D 


f-,INT,u 
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~INT,u 

F 

j?INT,p 


A  f 

e&Bh 


INT,u,e  fpINT,if  _ 
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e£Bh 


INT,<p,e 
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and 
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■  I  NT,  D 


■  INT,u 


-  INT,u 


- 1  NT,  ip 


A  (q(fINT’u’e)+qef1NT’u’e) 


e&Bh 


^ j;INT,u,e ^ 

J6e 


( cr )  dv 


- 1  NT,  ip 


A  ^^INT,p^+~efINT,p,e^ 
e£Bh 

^  I  NT,  u  i 


—  INT,D 

F  —  q 


~  INT,ip 


f/7VT,“=  a  (fNT’u’e) ,  fint,<?  =  a  (fINT^e) 


e£Bh 


e£Bh 


(119) 


- i  iV  i  ,  ~  . 

where  F  is  the  internal  force  and  moment  vector  associated  with  free  nodal  dofs  in  B, 

~  INT,D 

F  the  internal  force  and  moment  vector  associated  with  prescribed  nodal  dofs  in  B,  where 

~  INT,D  ~  INT,D  ~  /7VT,£> 

F  and  F  are  extracted  from  the  full  internal  force  and  moment  vector  F  ,  with 
superscript  (•) D  denoting  the  full  internal  force  and  moment  vector  associated  with  elements  in 
Bh,  qe  is  the  partitioning  coefficient  of  continuum  potential  energy  associated  with  element 
Be  C  Bh,  and  (•)  is  a  homogenization  operator  (to  be  defined  later).  A  simple  choice  is  a  volume 
fraction  qe  =  re.  (13)  considered  a  more  sophisticated  approach  using  an  atomic  bond  density 
function  solved  to  reproduce  a  minimum  potential  energy  state  for  homogeneous  deformation. 
The  analogy  here  for  particles  would  be  a  particle  contact  density  for  the  potential  energy  terms 
(internal  force  vectors).  This  will  be  considered  further  in  future  work.  For  now,  we  consider  a 
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volume  fraction  partitioning  through  qe,  and  a  simple  scaling  through  coefficients  q  (see  ID 

—  INT,D 

numerical  examples).  F  is  the  homogenized  internal  force  vector  associated  with 
prescribed  nodal  dofs  in  Bh,  which  has  no  contribution  if  q  =  0,  i.e.,  underlying  particle  contact 
forces  and  moments  provide  full  contribution  in  Bh.  FINT,D  is  the  internal  force  vector 
associated  with  free  nodal  dofs  in  the  pure  continuum  domain  Bh.  The  external  force  and 
moment  vectors  are  written  as 


j?EXT,D 


~EXT,D  pEXT,D 
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~  EXT,D 

where  F  is  the  external  body  force  and  couple  vector  associated  with  free  nodal  dofs  in  B  , 

~  EXT,D 

F  the  external  body  force  and  couple  vector  associated  with  prescribed  nodal  dofs  in  B  , 

~  EXT,D  ~  EXT,D  ~  EXT,D 

where  F  and  F  are  extracted  from  F  ,  the  total  external  body  force  and  couple 


—  EXT,D 

F  is  the  homogenized  external  body  force  and  couple  vector 


vector  calculated  in  Bh. 

associated  with  prescribed  nodal  dofs  in  Bh,  which  has  no  contribution  if  q  =  0,  i.e.,  underlying 
particle  body  forces  and  couples  provide  full  contribution  in  Bh.  F  '  ’  is  the  external  force 
and  couple  vector  associated  with  free  nodal  dofs  in  the  pure  continuum  FE  domain  Bh. 
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2.4.3  ID  Full  Overlap  Coupling 


Starting  with  equation  108,  and  referring  to  figure  7,  we  are  able  to  arrive  at  a  simplified  set  of 
nonlinear  equations  to  solve  for  a  fully  overlapped  particle-continuum  region. 


Figure  7.  Domain  of  full  overlap  coupling  between  a  1-D  string  of  twenty-one  glued,  Hertzian  nonlinear 


elastic,  spherical  DEs,  and  a  five-element  1-D  micropolar  linear  elastic  mesh. 

We  express  the  free  dofs  of  a  particle  a  in  terms  of  the  prescribed  nodal  dofs  of  an  element  e,  such 

that 
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where  £,(xa)  =  (2 /he)[xa  —  {x\  +  x%)/2]  ,  and  the  resulting  nonlinear  equations  to  solve  for  Q, 
with  D  =  BbqQ  and  FINT’b{D )  =  KD(BbqQ),  are 


(1  -  q)FINT’®(Q)  +  qBTbQKD(BbQQ )  =  FEXT,Q  (122) 


The  twenty-one  1-cm  diameter  DE  spheres,  and  overlaying  five-element  micropolar  elastic  mesh 
in  figure  7  are  used  to  demonstrate  the  overlap  coupling  procedure.  The  parameters  are  modified 
slightly  from  the  Timoshenko  beam  convergence  example:  FEXI  =  lOkN,  FEX1  =  lOkN, 

E  =  29GPa,  v  =  0.25,  k  =  0.1  /i,  and  7  =  m£2p,  £  =  1cm,  R  =  5mm,  B  =  50,  =  0.077, 

qv  =  5.5,  q9  =  1.0.  Notice  the  micropolar  couple  modulus  7  is  ~  3  *  400/0.05  =  24,  000  times 
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larger  than  used  for  the  results  reported  in  figure  4.  See  more  discussion  later.  The  different  q’s 
indicate  which  equation  they  influence  in  the  scaling  of  micropolar  elastic  FE  stiffness  (u  for 
axial,  v  for  transverse,  and  6  for  rotation).  The  axial  load  F^XT  is  first  applied  to  DE  particle  1 
(same  x  =  0  position  as  node  1  of  FE  mesh)  to  generate  a  transverse  stiffness  from  the 
Hertz-Mindlin  nonlinear  theory,  and  then  the  transverse  loading  FyXT  is  applied  while  the  axial 
load  is  held  fixed.  The  end  moment  at  x  =  0  is  zero.  It  is  cantilevered  at  x  =  L  (zero  axial  and 
transverse  displacement,  and  zero  rotation).  The  results  in  figure  8  demonstrate  how  the  DE  and 
FE  results  can  be  nearly  matched  by  scaling  the  micropolar  elastic  stiffness  with  the  g’s.  The 
rotations  of  the  FE  and  DE  do  not  match  closely,  and  there  is  no  coupling  between  DE  and  FE. 
This  example  demonstrates  the  scaling  through  q  of  the  micropolar  elastic  FE  stiffness.  Also, 
because  of  the  relatively  stiff  transverse  response  of  the  nonlinear  Hertzian  DE  spheres  in  contact 
after  1-cm  of  compression  (this  is  a  large  strain  compression  for  a  Hertzian  theory,  which  is  valid 
only  for  small  strains,  but  we  use  it  to  demonstrate  the  overlap  coupling),  the  value  of  the 
micropolar  couple  modulus  7  was  increased  by  a  factor  of  ~  3  *  400/0.05  =  24,  000,  thus  the 
large  transverse  displacement  observed  in  figure  4  versus  figure  8.  We  keep  the  parameters  of  the 
Hertzian  elastic  contact  model  constant,  and  vary  the  micropolar  elastic  model  parameters  to 
match  the  Hertzian  model.  This  micropolar  elastic  parameter  variation  and  energy  scaling 
through  the  q' s  take  the  place,  for  now,  of  the  homogenization  operator  (•)  discussed  in  the 
previous  section  on  the  overlap  coupling  method.  This  will  be  revisited  in  future  work.  Also,  in 
figure  8  we  use  10  times  the  axial  and  transverse  end  forces  F^XT  and  FyXT  when  compared  to 
figure  4. 

Using  the  same  micropolar  stiffness  scaling  coefficients  r/s,  and  equation  122,  and  with  energy 
factor  coefficients  qu  =  0.65,  qv  =  0.5,  qe  =  0.5,  we  get  the  results  in  figure  9.  Notice  that  the 
energy  is  factored  between  the  DE  particles  and  FE  micropolar  mesh.  The  DE  particle  rotations 
oscillate  about  the  FE  nodal  rotations.  Notice  that  the  particle  1  axial  force  ~  6.3kN  and  the 
particle  1  projected  from  rod  node  1  axial  force  ~  3.6kN  add  up  to  the  Ffx  1  =  lOkN  through 
equation  122.  Likewise  for  the  transverse  force.  This  is  the  effect  of  the  energy  factor 
coefficients  cf  s,  whereas  energy  partitioning  coefficient  q  has  no  effect  because  there  is  no  partial 
overlap  region  (i.e.,  Bh  =  0),  only  a  full  overlap  region  Bh  (see  figure  6). 

If  we  set  the  energy  factor  coefficients  qu  =  0 ,  qv  =  0,  q°  =  0,  we  get  the  results  in  figure  10. 
Notice  that  the  micropolar  FE  mesh  contributes  no  energy  to  the  system.  The  particle  DE 
simulation  provides  all  the  energy,  and  in  the  rotations,  they  oscillate. 
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AXIAL  FORCE  (kN)  AXIAL  DISPLACEMENT  (cm) 


Figure  8.  Full  overlap  results  with  scaling  coefficients  q’s  on  micropolar  stiffness,  but  no  coupling. 
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AXIAL  DISPLACEMENT  (cm) 


Figure  9.  Full  overlap  coupling  results  with  energy  factor  coefficients  qu  =  0.65,  qv  =  0.5,  qe  =  0.5. 
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Figure  10.  Full  overlap  coupling  results  with  energy  factor  coefficients  qu  =  0,  qv  =  0,  q°  =  0. 
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2.4.4  ID  Partial  Overlap  Coupling  with  Partial  Overlay  ID  Micropolar  FE 


Referring  to  figure  1 1 ,  NqD  =  0  because  there  is  no  overlap  between  free  particle  dofs  and  free 
continuum  dofs.  We  need  the  additional  interpolation  matrices  for  particle  a 
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where  we  can  then  relate  prescribed  dof  to  free  dof  through  the  projection  operators  as 
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where  we  solve  the  coupled  system  of  nonlinear  balance  equations  as  (simplified  from  equations 
108  and  109) 


Rq(Q,D)  =  Fint’q(Q)  +  B^qFint’Q{Q)  +  B^qFint’S{D)  -  Fext’q  =  0 


DQ 
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Figure  11.  Domain  of  partial  overlap  coupling  between  a  ID  string  of  1 1  glued,  Flertzian  nonlinear  elastic, 
spherical  DEs,  and  a  4  element  ID  micropolar  linear  elastic  mesh. 

For  volume  average  energy  partitioning  in  the  overlap  region  Bh, 

q  =  ( heAe  —  2.5(4/3)7ri?3)/ (heAe)  =  0.583  because  there  are  2.5  particle  volumes  in  Bh.  We 
illustrate  the  performance  of  the  overlap  coupling  algorithm  in  figures  12  and  13.  We  can  see  that 
with  scaling  of  micropolar  elasticity  through  the  q  coefficients,  we  can  achieve  a  homogeneous 
axial  displacement  gradient  across  the  overlap  coupling  region  Bh,  whereas  the  transverse 
component  cannot  be  made  homogeneous.  This  can  be  observed  because  of  the  ratcheting  of  the 
DE  particle  rotations,  which  the  micropolar  continuum  FE  can  only  represent  if  the  element 
length  he  is  chosen  to  be  one  DE  particle  diameter,  which  defeats  the  purpose  of  the  overall 
overlap  coupling  strategy.  We  believe  that  with  more  particles  in  2D  and  3D,  as  illustrated  in 
figure  6,  the  particle  displacements  and  rotations  will  be  smoothed  and,  thus,  the  overlap  coupling 
should  work  more  effectively.  In  a  sense,  this  ID  example  assumes  “too  discrete”  a  ID  DE 
particle  string,  that  exhibits  an  oscillation/ratcheting  behavior  upon  transverse  shear  loading.  The 
axial  component  is  handled  without  trouble. 

We  set  the  micropolar  scaling  coefficients  q u  =  0.0195,  q°  =  0.8,  qe  =  0.85  in 
Bh  =  Bh  U  Bh  U  Bh  for  ID  micropolar  FE  stiffness;  energy  factor  coefficients  qu  =  0,  qv  =  0, 
qe  =  0  in  Bh  U  Bh ;  and  energy  partitioning  coefficient  q  =  0.583  in  Bh.  The  “not  scaled”  plots 
indicate  that  qu  —  1,  q"  —  1,  q8  —  1.  The  results  are  shown  in  figure  12. 

If  we  set  the  micropolar  scaling  coefficients  the  same  q™  =  0.0195,  q"  =  0.8,  q8  =  0.85,  and 
energy  factor  coefficients  qu  —  0,  qv  —  0,  qe  =  0.9,  then  we  have  the  result  in  figure  13.  Notice 
how  factoring  some  of  the  DE  rotational  energy  to  the  ID  micropolar  FE  mesh  reduces  the 
oscillations  of  the  DE  particles  in  the  transition/overlap  region. 
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AXIAL  DISPLACEMENT  (cm) 


Figure  12.  Partial  overlap  coupling  results  for  energy  factor  coefficients  qu  =  0,  qv  =  0,  q°  =  0. 
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Figure  13.  Partial  overlap  coupling  results  for  energy  factor  coefficients  qu  =  0,  qv  =  0,  qe  =  0.9. 
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3.  Summary 


3.1  Results 

The  details  of  a  1-D  overlap  coupling  between  a  micropolar  linear  isotropic  elastic  finite  element 
(FE)  model  and  a  ID  string  of  Hertzian  (nonlinear)  elastic-at-contact  discrete  element  (DE) 
spheres  were  presented.  Numerical  examples  demonstrated  various  overlapping  domains:  (1) 
full  overlap  coupling  with  fully-informed  upscaled  micropolar  FE  response  from  the  underlying 
DE  response,  and  (2)  partial  overlap  coupling  with  partially-informed  upscaled  micropolar  FE 
response  (along  partial  length  of  domain)  from  the  underlying  DE  response  (also  along  partial 
length  of  domain). 

3.2  Conclusions 

The  simple  ID  problem  presented  in  this  report  provides  an  accessible  model  problem  through 
which  to  better  understand  how  such  coupling  strategies  should  work  for  overlap  coupling  of 
underlying  DNS  particulate  models  (in  this  case,  DE),  and  overlying  generalized  continuum 
models  (in  this  case,  micropolar  elasticity).  What  makes  this  coupling  strategy  different  than 
those  for  atomistic  continuum  coupling  methods  (13,  19)  is  the  rotational  degrees  of  freedom  of 
the  DE  model,  the  open  window  BDE  region,  and  the  additional  degrees  of  freedom  inherent  in 
the  generalized  continuum  models — in  this  case  a  micropolar  continuum  (a  subset  of 
micromorphic  continuum). 

3.3  Future  Work 

Future  work  will  first  involve  extending  the  ID  formulation  and  implementation  to  include  inertia 
terms  to  study  wave  propagation  (axial,  transverse,  and  rotational)  along  the  ID  domain. 
Extending  these  concepts  for  the  envisioned  full  3D  overlap  coupling  problem  and  micromorphic 
continuum  as  outlined  in  section  1.1  will  come  next. 
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